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ABSTRACT 

When clusters of galaxies are viewed in projection, one cannot avoid picking up a fraction of foreground/background interlopers, that 
f \ ' lie within the virial cone, but outside the virial sphere. Structural and kinematic deprojection equations are known for the academic 

case of a static Universe, but not for the real case of an expanding Universe, where the Hubble flow (HF) stretches the line-of-sight 
, distribution of velocities. Using 93 mock relaxed clusters, built from the dark matter (DM) particles of a hydrodynamical cosmological 

i' simulation, we quantify the distribution of interlopers in projected phase space (PPS), as well as the biases in the radial and kinematical 

structure of clusters produced by the HF. The stacked mock clusters are well fit by an m=5 Einasto DM density profile (but only out to 
1 .5 virial radii), with velocity anisotropy (VA) close to the Mamon-Lokas model with characteristic radius equal to that of density slope 
-2. The surface density of interlopers is nearly flat out to the virial radius, while their velocity distribution shows a dominant gaussian 
cluster-outskirts component and a flat field component. This distribution of interlopers in PPS is nearly universal, showing only small 
trends with cluster mass, and is quantified. A local k=2.1 sigma velocity cut is found to return the line-of-sight velocity dispersion 
profile (LOSVDP) expected from the NFW density and VA profiles measured in three dimensions. The HF causes a shallower outer 
, LOSVDP that cannot be well matched by the Einasto model for any value of k. After this velocity cut, which removes 1 interloper out 

of 6, interlopers still account for 23+1% of all DM particles with projected radii within the virial radius (surprisingly very similar to 
' the observed fraction of cluster galaxies lying off the Red Sequence) and over 60% between 0.8 and 1 virial radius. The HF causes the 

| best-fit projected NFW or m=5 Einasto model to the stacked cluster to underestimate the true concentration measured in 3D by 6+6% 

, (16+7%) after (before) the velocity cut. These biases in concentration are reduced by over a factor two once a constant background is 

■ included in the fit. The VA profile recovered from the measured LOSVDP by assuming the correct mass profile recovers fairly well 

the VA measured in 3D, with a slight, marginally significant, bias towards more radial orbits in the outer regions. These small biases 
in the concentration and VA of the galaxy system are overshadowed by important cluster-to-cluster fluctuations caused by cosmic 
variance and by the strong inefficiency caused by the limited numbers of observed galaxies in clusters. An appendix provides an 
analytical approximation to the surface density, projected mass and tangential shear profiles of the Einasto model. Another derives the 
expressions for the surface density and mass profiles of the NFW model projected on the sphere (for future kinematic modeling). 



Key words. Galaxies: clusters: general — Cosmology: miscellaneous — (Cosmology): dark matter — Galaxies: halos — 
Gravitational lensing: weak — Methods: numerical 

1 . Introduction where S and v are the projected and space number densities, re- 
spectively, while R and r are the projected and space radial dis- 

The galaxy number density profiles of groups and clusters of tances (hereafter, radii), respectively. Equation (1) can be depro- 

galaxies falls off slowly enough at large radii that material be- j ec ted through Abel inversion 1 to yield 

yond the virial radius (within which these structures are thought „ 

to be in dynamical equilibrium) contribute non-negligibly to the v ( r ) _ __ I dR (2) 

projected view of cluster, i.e. to the radial profiles of surface x J r \Ir 2 - r 2 

density, line-of-sight velocity dispersion and higher velocity mo- The jection t0 infinit is licit in equations (1) and (2 ). 

However, the Hubble expansion complicates the picture, as 
In principle, this contamination of observables by interlop- the Hubble flow moves background (foreground) objects to high 
ers, defined here as particles that lie within the virial cone but positive (negative) line-of-sight velocities. This is illustrated in 
outside the virial sphere, is not a problem, since we know how Figure 1 which shows how the line-of-sight velocity vs. real- 
to express deprojection equations when interlopers extend to in- space distance relation is affected by the Hubble flow. The line- 
finity along the line-of-sight. Consider the projection equation of-sight distances are computed as the segment length QP in 

+ra ^ 1 Alternatively, the projection equation (1) corresponds to a con- 

_ f" y(r)ds — 2 f v(r) r ^ r (1) v °l ut i° n an d can therefore be deprojected with Fourier methods (see 

J m Jr VV 2 — R 2 ' Discussion in Mamon & Boue 2010 and references therein). 
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Fig. 1. Line-of-sight velocity as a function of real-space line-of-sight 
distance (see Fig. 2) for particles inside the virial cone obtained by 
stacking 93 cluster-mass halos in the cosmological simulation described 
in Sect. 2 without (top) and with (bottom) the Hubble flow (1 particle in 
5 is shown for clarity). The red dashed horizontal lines roughly indicate 
the effects of a radius-independent 3 cr clipping. Note that the velocity- 
distance relation without Hubble flow shown here is not entirely real- 
istic, because the simulation was run in the context of an expanding 
Universe (and cannot be run in a static Universe, for lack of knowledge 
of realistic initial conditions), but should be accurate enough to illus- 
trate our point. 



Figure 2. Now, clipping the velocity differences to, say, k — 3 
times the cluster velocity dispersion (averaged over a circular 
aperture, hereafter aperture velocity dispersion), o~ v , gets rid of 
all the distant interlopers. More precisely, the radius, r max , where 




Fig. 2. Representation of the virial cone with halo particles inside the 
inscribed virial sphere and interlopers outside. Also shown is our defi- 
nition of projected radius (CQ) and line-of-sight distance (OP and QP 
respectively in the observer and halo reference frames) for a random 
point P. For illustrative purposes, the distance to the cone is taken to be 
very small, so that the cone opening angle is much larger than in reality. 



the Hubble flow matches k <x„ is found by solving Hq r max = kct u 
(where Ho is the Hubble constant) yielding 




(3) 



where r v is the virial radius where the mean density is A times 
the critical density of the Universe, p c = 2>H^I(%nG) (where G 
is the gravitational constant), and where v v = \[AJ2Hq r v is the 
circular velocity at the virial radius. Clusters are thought to have 
density profiles consistent with the Navarro, Frenk, & White pro- 
file (1996, hereafter NFW) profiles, 



v(r) 
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(4) 



where r_2 is the radius of density slope -2, in number (Lin, 
Mohr, & Stanford 2004), luminosity (Lokas & Mamon 2003) 
and mass (Lokas & Mamon; Biviano & Girardi 2003; Katgert, 
Biviano, & Mazure 2004), with a concentration, c = r v /r_2, of 3 
to 5. For isotropic NFW models, the aperture velocity dispersion 
is (Appendix A of Mauduit & Mamon 2007) cr„ = rj v v , where 
77 ^ 0.62 (weakly dependent on concentration), and equation (3) 
then becomes 



T JSL * 13.2 (*) J± 
r v \3) V 100 



(5) 



Equation (5) indicates that a 3-sigma clipping will remove all 
material beyond 13 (A = 100) or 19 (A = 200) virial radii. 2 
So, in the deprojection equation (2), the upper integration limit 
must be set to this value of r max . Although this effective cutoff 
in line-of-sight distances is quite far removed from the cluster, 
it is not clear whether there may still be a measurable bias in 
the concentration of clusters that one measures by comparing 
the surface density distribution of galaxies in clusters with NFW 



2 With our chosen cosmology, the overdensity at the virial radius is 
A a 100, but many authors prefer to work with A = 200, and we will 
do so too. 
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models projected out to infinity. Moreover, it is not clear how 
accurate are such measures given the finite number of galaxies 
observed within clusters. 

Finally, it is not clear how the stretching of the velocities af- 
fects the kinematic analyses of clusters, especially in the case 
of nearby clusters where the opening angle of the cone is non- 
negligible, leading to an asymmetry between the foreground and 
background absolute velocity distributions. For example, is the 
anisotropy of the 3D velocity distribution (hereafter velocity 
anisotropy or simply anisotropy) 



i-I 

2 



(6) 



or equivalently 3 



W) = [i -ffi] 



-1/2 



2 (iftr)) 



1/2 



(7) 



affected by the Hubble flow? One can add interlopers beyond 
the virial radius as a separate component to the kinematical 
modeling (van der Marel et al. 2000; Wojtak et al. 2007). 
Unfortunately, we have no knowledge of the distribution of in- 
terlopers in projected phase space (projected distance to the halo 
center and line-of-sight velocity). 

This paper provides the distribution of interlopers in pro- 
jected phase space (projected distance to the halo center and line- 
of-sight velocity) as measured on nearly 100 stacked halos from 
a well-resolved cosmological simulation. We additionally mea- 
sure the bias in the measured surface density and line-of-sight 
velocity dispersion and kurtosis profiles compared to those ob- 
tained in a Universe with no Hubble flow, and estimate how this 
bias affects the recovered concentration and velocity anisotropy 
of the cluster. In this paper, we use interchangeably the terms 
'clusters' and 'halos'. 

We present in Sect. 2 the cosmological simulations we use 
and how the individual halos were built. In Sect. 3 we explain 
how we stack these halos. In Sect. 4, we present the statistics on 
the halo members and the interlopers in projected phase space. 
Then in Sect. 5, we explain how we remove the outer interlop- 
ers and show analogous statistics on the cleaned stacked halo in 
Sect. 6. We proceed in Sect. 7 to measure the biases induced by 
the Hubble flow and the imperfect interloper removal on the es- 
timated concentration parameter and anisotropy profile. We dis- 
cuss our results in Sect. 8. 



and entropy conservation, radiative cooling, a uniform time- 
dependent UV background (Haardt & Madau 1996), the self- 
regulated hybrid multi-phase model for star formation (Springel 
& Hernquist 2003), and a phenomenological model for galactic 
winds powered by Type-II supernovae. 

Dark matter halos were identified by Borgani et al. at red- 
shift z — by applying a standard Friends-of-friends (FoF) 
analysis to the dark matter particle set, with linking length 0.15 
times the mean inter-particle distance. After the FoF identifica- 
tion, the center of the halo was set to the position of its most 
bound particle. A spherical overdensity criterion was then ap- 
plied to determine the virial radius, r v = r^xi of each halo. In 
this manner, 117 halos were identified within the simulated vol- 
ume, among which 105 form a complete subsample with virial 
mass M200 larger than 10 14 /i _1 M o , thus representing a sample 
of mock galaxy clusters. Their mean and maximum masses are 
respectively 2.0 x 10 14 /z -1 M o and 1.1 x 1O 15 ^ _1 M . 

To save computing time, we worked on a random subsample 
of roughly 2 million particles among the 480 3 . Although the sim- 
ulation also produced galaxies, we chose to use the dark matter 
particles as tracers of the galaxy distribution for two reasons: 1) 
Simulated galaxy properties in cosmological simulations depend 
on details of the baryon physics implemented in the code, and 
can show some mismatch with observed properties (e.g. Saro 
et al. 2006); 2) only a handful of simulated clusters had over 
50 galaxies (Saro et al.), so we would have strongly suffered 
from small-number statistics. There is some debate on whether 
the velocity distribution of galaxies is biased relative to the dark 
matter. On one hand, the galaxy velocity distribution, although 
close to the dark matter one, shows a preference for lower ve- 
locities (Biviano et al. 2006), perhaps as a consequence of dy- 
namical friction. On the other hand, the velocity distribution of 
subhalos, selected with a minimum mass before entering their 
parent cluster-mass halos, is similar to that of the dark matter 
(Faltenbacher & Diemand 2006). 

We visually inspected each of the 105 clusters in redshift 
space along three orthogonal viewing axes, and removed 12 clus- 
ters that appeared, within r2oo> to be composed of two or three 
sub-clusters of similar mass (where the secondary had at least 
40% of the mass of the primary). Most observers would omit 
such clusters when analyzing their radial structure or internal 
kinematics. This leaves us with 93 final mock clusters. 4 The 
median values (interquartile uncertainties) of their virial radii, 
virial masses, virial circular velocities, and velocity dispersions 
(within their virial spheres) are respectively 864+81 h~ l kpc, 
1.50±0.44 \0 u h- l M Q , 865+81 kms- 1 , and 584+60^8-'. 



2. Data from cosmological TY-body simulations 

The halos analyzed in this paper were extracted by Borgani et al. 
(2004) from their large cosmological hydrodynamical simula- 
tion performed using the parallel Tree+SPH GADGET-2 code 
(Springel 2005). The simulation assumes a cosmological model 
with present day parameters Q m = 0.3, Qa = 0.7, Q.^ = 0.039, 
h = // /(100 km s-'Mpc- 1 ) = 0.7, and cr 8 = 0.8. The box 
size is L — 192 h~ l Mpc. The simulation used 480 3 dark mat- 
ter particles and (initially) as many gas particles, for a dark mat- 
ter particle mass of 4.62 x 10 9 lr l M & . The softening length was 
set to 22.5 comoving kpc until z = 2 and fixed afterwards 
(i.e., 7.5 h~ x kpc). The simulation code includes explicit energy 



3 (3(r) enters the Jeans equation of local hydrostatic equilibrium, 
while !Mf) is a more physical definition of velocity anisotropy. 



3. Stacking the virial cones 

For each cluster, we projected the coordinates along the virial 
cone (circumscribing the virial sphere, see Fig. 2), as follows. 
We first renormalized the 6 coordinates of phase space of the 
entire simulation box to be relative to the cluster. So the cluster 
most bound particle should be at the origin and its mean peculiar 
(bulk) velocity should be zero. 5 To take into account the peri- 
odic boundaries of the simulation box, we added or subtracted 
a box length to those particles situated at over a half -box length 



4 Including all 105 halos makes virtually no difference for the results 
of this article. 

5 Observers usually adopt the position of the brightest cluster galaxy 
as the center, and this corresponds to the most bound galaxy, so they 
should not suffer from important centering errors, although admittedly 
some clusters like Coma have two brightest galaxies. 
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from the cluster center. In this fashion, each cluster now effec- 
tively sits at the center of the simulation box. We then placed an 
observer at coordinates (-D, 0, 0), (0, —D, 0), or (0, 0, -D), with 
< D < L/2. We present here the results for D - 90 /T 1 Mpc, 
corresponding to a typical distance of observed clusters in the 
local Universe. At this distance, the median virial angular radius 
of our 93 clusters is 33+3 arcmin. We do not expect that the re- 
sults of this paper should depend on the adopted value of D. We 
assume that the observer's peculiar velocity is equal to the clus- 
ter's bulk velocity, so that the observer's velocity is zero in the 
renormalized coordinate system. 6 

We then measured, for each cluster and for each of these 
three observers, the coordinates of all 2 million particles in both 
the observer frame and the cluster frame. Given the distance r D 
of the particle to the observer and the projected coordinate of the 
particle in cylindrical coordinates R c , we determined the pro- 
jected distances in the observer frame (measured in a plane per- 
pendicular to the line-of-sight passing through the cluster cen- 
ter) as R = DRJr (see Fig 2, where D = OC, r = OP, 
R c = PP', and R a = CQ). This projection ensures that parti- 
cles along the surface of the virial cone have R — r v . We were 
then able to select all particles within a cone circumscribing 
the sphere of radius r v , where we chose r w = r2oo> as we ll as 
r v = 1.35 r2oo — Hoo- In practice, we extracted data from a wider 
cone, circumscribing the sphere of radius 3 f2oo- 

We next added the Hubble flow (for both the observer and 
cluster frames), using Hq = 100 km s _1 Mpc -1 ? We then limited 
in depth to line-of-sight velocities within 4 v v from the cluster. 
As mentioned in Sect. 1, the Hubble flow effectively limits the 
depth of the cones to a half-length of k y/A/2 = 40 virial radii 
(see eq. [3]), where the cut in velocity space in units of virial 
velocity is k — 4 and the virial overdensity is A = 200. 8 Our 
results should not depend much on the distance of our observer 
(90/0.864 = 104 virial radii), as our first cut at +4 virial veloci- 
ties limits the line of sight to 40% of the observer's distance. 

We finally normalized the particle relative positions and ve- 
locities to the virial radius, r v and circular velocity at the virial 
radius, v v , respectively, and finally computed the projections of 
the velocities in spherical coordinates (to later measure the 3D 
radial profiles of density and velocity anisotropy). For clarity, we 
will sometimes use the notation f-200 for the virial radius and D200 
for the virial velocity. 

In the end, we have roughly 84 thousand particles for each 
of the three cartesian stacked virial cones, which we then stack 
together into our global stacked virial cone, with a grand total of 
roughly a quarter-million particles, among which nearly three- 
quarters lie within the virial radius. Hence, roughly l/30th of 
all the particles in the simulation box lie within the virial radius, 
r2oo, of our 93 clusters. Note that with our stacking method, some 
of the particles inside virial cones but outside virial radii can 
end up being selected in more than one of the three cartesian 
stacked virial cones. However, this fraction is small (27% of the 
interlopers, i.e. less than 8% of all particles in virial cones), so 
the three cartesian stacked virial cones are virtually independent 



6 In our simulation, the one-dimensional cluster bulk velocity disper- 
sion is 228kms~', so given our adopted distance of D = 90 h Mpc, 
the typical cluster bulk velocity is only 228/9000 = 2.5%. Therefore, 
our neglect of the observer's peculiar motion relative to the cluster is an 
adequate assumption. 

7 All our positions and cluster virial radii were expressed in /T 1 kpc; 
the choice of Hq does not matter as long as we normalize to the virial 
quantities r v and u v , which we computed with the same value of Hq. 

8 Throughout this paper, we use x to express quantity x in virial units. 



(except that their halos are common), hence can be stacked into 
the global virial cone. 
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Fig. 3. Projected phase space diagram of stacked virial cone (built from 
3 x 93 halos). Only 1 particle in 5 is shown for clarity. The red (lighter) 
and blue (darker) points refer to the particles within and outside the 
virial sphere, respectively. The green curves illustrate the ±2.1 <T\ 0! ,(R) 
velocity cut (from eqs. [21] and [22]) for the c = C200 = 4 NFW model 
with r a =r_2 ML anisotropy (eq. [20]). 

Figure 3 shows the projected phase space (line-of-sight ve- 
locity, Vios, versus projected radius, R) distribution of particles of 
the stacked cluster, highlighting the halo particles (lying within 
the virial sphere) and interlopers (lying in the virial cone but out- 
side of the virial sphere). Note that the halo particles are confined 
to fairly small velocities (the largest absolute halo particle veloc- 
ity is 2.7 v v ). One notices an excess of positive velocity outliers 
in comparison with negative velocity ones, as expected from the 
conical projection used here. 

4. Interloper statistics before the velocity cut 

We now measure the distribution of interlopers in projected 
phase space and study its dependence on halo mass. 

4.1. Global statistics 

Figure 4 displays contours of the density in (R, \i>i os \) projected 
phase space. Note that the projected phase space density of 
Figure 4 is proportional to d 2 A7(/?d/?d|f;i os |), hence the differ- 
ent shapes than seen in Figure 3. The interlopers have a very 
different projected phase space density than the halo particles. 
In particular, their horizontal contours mean that the interloper 
projected phase space density is fairly independent of projected 
radius. Moreover, the interloper contours do not extend beyond 
uios = 2 u v , except for a few islands (caused by cosmic variance), 
indicating that the velocity distribution of interlopers is close to 
flat at large velocities (the islands thus represent small, probably 
statistical, fluctuations in a flat background). 
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Fig. 4. Contours of projected phase space density of stacked virial cone 
— in units of d 2 N/(dR/r2oo d|i'iosl/f2Qo)/(2?r/f/r2oo) — of halo (dashed 
red contours) and interloper (solid, blue) particles. The halo contours 
are equally spaced in log-density, increasing by a factor 2.44 from 0.085 
(upper right) to 3370 (lower left). The interloper contours are taken from 
the same set as the halo contours, but limited from the 5th highest level 
halo contour (72 virial units, at the lower left) to the 8th highest level 
halo contour (4.0 virial units, for all contours above the 3rd nearly hor- 
izontal one). The green curve shows the 2.7 cr los (R) (from eqs. [21] and 
[22]) velocity cut for the c-4 NFW model with r a =r_2 ML anisotropy 
(eq. [20]). 




These issues can be looked in more detail through slices of 
the projected phase space density in velocity and radial space. 
Figure 5 shows how the projected phase space density varies 
with projected radius in different wide velocity bins. The halo 
particles display a negative gradient, i.e. a decreasing surface 
number density profile, as expected. However, one immediately 
notices that in all line-of-sight velocity bins, the density of inter- 
lopers in projected phase space is roughly independent of pro- 
jected radius. In other words, interlopers have a nearly flat sur- 
face density profile. 

For low velocities, one can notice a small rise of the inter- 
loper surface density at high projected radii. This small rise is 
a geometric effect: the line-of-sight distance between the virial 
sphere and a sphere of k > 1 virial radii isle - I virial radii at R = 
but VA: 2 - 1 virial radii at R = r y , which is y/(k + l)/(k- 1) 
times greater. We will return to this rise in Sect 4.2. 

Figure 6 displays the distribution of line-of-sight velocities 
of the halo and interloping particles. The interloper distribution 
shows a flat component that dominates at large velocities and 
a gaussian-like component. Figure 6 confirms, once more (see 
Figs. 4 and especially 5) that the density of interlopers in pro- 
jected phase space is fairly independent of radius. 

The total surface density of particles in velocity space 
shows an inflection point at about 2 v v (bottom plot of Fig. 6). 
Kinematical modelers attempt to throw out the high-velocity in- 
terlopers by identifying this gap by eye (Kent & Gunn 1982; 
Lokas & Mamon 2003) or automatically (Fadda et al. 1996) or 
by rejecting 3 cr outliers, either using a global criterion (Yahil 



Fig. 5. Phase space density of stacked virial cone as a function of pro- 
jected radius in bins of absolute line-of-sight velocities (marked on the 
upper-right of each plot). Dashed (red) and dotted (blue) histograms 
represent the halo particles (r < r Y ) and the interlopers (r > r v ), respec- 
tively, while the solid histograms (artificially moved up by 0.06 dex for 
clarity) represent the full set of particles. There are no halo particles at 
v > 3 v v (top plot). 



& Vidal 1977) or a local one (e.g. Lokas et al. 2006; Wojtak & 
Lokas 2010). Interestingly, the 3 <x„ criterion was first motivated 
on statistical grounds, but the v 2.0 v w inflection point one sees 
in the plots of Figure 6 happens to correspond to 3 cr„. In other 
words, the 3 <x„ criterion is not only a consequence of statistics, 
but also motivated by the combination of cluster dynamics and 
cosmology. 

While visual attempts to separate interlopers from halo par- 
ticles in projected phase space look for gaps in the line-of-sight 
velocity distribution, the different panels of Figure 6 indicate 
that, on average, one should not expect such gaps in the pro- 
jected phase space diagram of stacked clusters, as the number 
of interlopers also decreases with velocity to reach a plateau at 
about 2 iv 
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Fig. 6. Phase space density of stacked virial cone as a function of line- 
of-sight velocity in different radial bins (marked on the lower-left of 
each plot). Dashed ( red ) and dotted ( blue ) histograms represent the halo 
particles {r < r v ) and the interlopers (r > r v ), respectively, while the 
solid histograms (artificially moved up by 0.06 dex for clarity) represent 
the full set of particles. The green vertical lines indicate the 2.7 cr los {R) 
(from eqs. [21] and [22]) velocity cut for the c=4 NFW model with 
r il =r_2 ML anisotropy (eq. [20]). 



4.2. Universality 

The three panels of Figure 7 illustrate the universality of the line- 
of-sight velocity distribution of interlopers, in terms of projected 
radius and halo mass. Figure 7a confirms that the projected phase 
space density of interlopers depends little on projected radius. 
This will be quantified later in this sub-section. 

As mentioned above, to first order, the density of interlopers 
in projected phase space has a constant component and a quasi- 
gaussian component, which we write 



g(R, Ksl) - A exp 



+ B . 



(8) 



Maximum likelihood estimation (MLE, see appendix C) yields 
cri = 0.576 ± 0.003 , 
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Fig. 7. Phase space density of interlopers as a function of line-of-sight 
velocity, a {top): Dependence on projected radius: R/Pm = 0-0.2, 0.2- 
0.4, 0.4-0.6, 0.6-0.8, and 0.8-1 {red dotted, green short-dashed, blue 
long-dashed, magenta dot-long-dashed and cyan dot-long-dashed bro- 
ken lines, respectively), b {middle): Dependence on 3D radial distance: 
1 < r/rjoo < 8 {blue dash-dotted curve) and r/r2ao > 8 {red dashed 
curve), c {bottom): Dependence on halo mass: high (/7M200 > 1.87 x 
10 14 M Q , red long dashed curve) and low {h M 2 oo < 1.87x 10 14 M e , blue 
dash-dotted curve) mass halos. Also shown in all three plots is the MLE 
(eqs. [8] and [10]) for the gaussian {dotted curve) and field {horizontal 
dashed line) components, as well as the predicted total interlopers (the 
sum of these two components, solid curve). 



A = 0.1164 + 0.0006, 
B = 0.0075 + 0.0001 , 



(9) 



where cr\ is in units of v v , while g, A and B are in units of 
N v r~ 2 v~ , where N v is the number of particles within the virial 
sphere. Figures 7a and 7c show that this gaussian+constant 
model is an excellent fit to the distribution of interloper line- 
of-sight velocities. 

The origin of these two components is clarified in Figure 7b: 
cutting the halo interlopers in two subsamples at different 3D 
distances from the halo, one finds that the flat component cor- 
responds roughly to the interlopers beyond 8 r%x>, while the 
quasi-gaussian component corresponds to the closer interlopers 
O200 < r < 8 r 2 oo). 

Figure 8 shows the radial dependence of the parameters of 
the interloper phase space distribution. We also show the mean 
surface densities of interlopers measured on the stacked halo. 
The normalization of the gaussian component and the surface 
density both increase slowly at small projected radii, but sharply 
near the virial radius. A good fit for the normalization and the 
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Fig. 8. Variations of MLE parameters of equation (8) and measured 
interloper surface density with projected radius (in units of v v for en, 
N y r~ 2 v~ l for A and N v r~ 2 for £;). The blue and green curves show the 
fits of equations (10) and (11), respectively, while the red horizontal 
dashed line shows the mean of 50 B. The magenta solid curve shows 
the prediction for the interloper surface density from equations (12), 
(10) and (11), respectively, with B = 0.0075, while the magenta dashed 
curve shows the interloper mean surface density predicted (eq. [13]) 
from a c=4 NFW model. Small radial bins were chosen to capture the 
rise of Ej and A near the virial radius. Errors are from the likelihood 
ratios in the MLE fit (A, o~\, and B) or Poisson (Si). 



standard deviation of the gaussian component is provided by 

A(R) dex [- 1. 06 1 +0.364 X 2 - 0.580 X 4 + 0.533 X 6 ] , (10) 
o-i(R)^ 0.612- 0.0653 X 2 , (11) 

where X = R/R200 and dex x =10*. 

The solid magenta curve of Figure 8 indicates that the inter- 
loper surface density profile is well recovered from A, cr\ and B 
by integrating over the model velocity distribution (eq. [8]) from 
to kv v : 



£;(/?) = V^/2A(^)cr i (/?)erf 



ri(R) V2 



+ kB . 



(12) 



Alternatively, the surface density profile of interlopers is also 
well recovered (dashed magenta curve in Fig. 8) by the predic- 
tion from an NFW model, i.e. as the difference between the stan- 
dard surface density integrated to infinity and the surface density 
limited to the virial sphere: 



Ej(/?) - S NFW (^) 



(13) 



where the expression for S^p W (/?,7?20o) is provided in equa- 
tion (B.4). According to equation (13), the rise at radii close 
to the virial radius is almost independent of the concentration 
(with relative differences of less than 4% at all projected radii). 
Hence, one can adopt a unique empirical model for A. These two 
estimates of E;(R) are in very good agreement, except that the 
difference of the cylindrical and spherical NFW surface density 



profiles (eq. [13]) rises somewhat faster than our fit at projected 
radii very close to the virial radius. 

The bottom panel of Figure 7 shows that the density of in- 
terlopers in projected phase space as a function of line-of-sight 
velocity remains the same for low and high mass clusters (where 
we took the dividing line at the median cluster virial mass of 
/1M200 = 1.87 x 1O 14 M ). The interlopers of high mass clusters 
have a cluster-outskirts component whose density in projected 
phase space is roughly 15% lower than that of the low-mass clus- 
ters. 




10 14 



2xl0 14 
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Fig. 9. Variations of MLE parameters of equation (8) with halo mass (in 



units of v y for o~\, N v r v for £j and 7V V r~ 



for A and B). Errors are 



from 100 bootstraps on the 93 halos. The symbols for 50 B are displaced 
by +0.01 dex for clarity. 



Figure 9 provides a closer look at the variation with halo 
mass of the parameters of equation (8). The best fitting logarith- 
mic slopes, obtained by least-squares fits to the points shown 
in Fig. 9 are -0.042+0.030, -0.037+0.069, 0.155+0.225 and 
-0.044+0.104 for o- u A, B, and respectively. The 90% con- 
fidence lower limits on the slopes dlnX/dlnM are thus -0.08, 
-0.13, -0.17 and -0.18 for X = a u A, B, and S;, respectively, 
while the 90% upper limits are less than 0. 1 except for B (where 
is it 0.5). These shallow limits to the logarithmic slopes illus- 
trate the near-universality of the distribution of interlopers in 
projected phase space. 

However, the universality of the distribution of interlopers 
in projected phase space may hide an important level of cosmic 
variance. Performing MLE for parameters en, A, and B, for each 
of the 93 clusters, each viewed in turn along each of three or- 
thogonal viewing axes, we find standard deviations of 



0-(logo-i) = O.ll, cr(logA) = 0.23, o-(logB) = 0.40. (14) 

So, while the dispersion en of the gaussian component of the in- 
terloper velocity distribution is fairly constant (29% typical vari- 
ations) from one cluster to the next, there is more scatter in the 
normalization A of the gaussian component (factor 1 .7 typical 
variations) and a large scatter in the flat component B (factor 2.5 
typical variations). 
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5. Interloper removal 

We remove the interlopers of the stacked cluster proceeding 
along similar lines as Lokas et al. (2006) (see also Wojtak et al. 
2007), by clipping the velocities beyond k times the local line- 
of-sight velocity dispersion. We assume that our stacked cluster 
has an NFW profile (eq. [4]), or, alternatively, an Einasto (1965) 
profile: 



y(r) = 



(2m) 



3m 



m y(3m, 2m) 



exp 



—2m — 



l/m 



M(r_ 2 ) 



(15) 



where y(a,x) is the incomplete gamma function. The density 
model of equation (15) fits the density profiles of ACDM halos 
even better than the NFW model (as first discovered by Navarro 
et al. 2004), at the expense of an additional parameter, m. 

The velocity cut requires an estimate of the line-of-sight ve- 
locity dispersion profile of the halo. One could measure this in 
bins of projected radius, iteratively rejecting the outliers. This 
gives a profile that shows important radial fluctuations, and we 
would need to either smooth the profile or fit a smooth analytical 
function to it. 

Instead, we choose to predict the line-of-sight velocity dis- 
persion profile given the typical density and velocity anisotropy 
profiles of halos. The line-of-sight velocity dispersion profile can 
be written (Mamon & Lokas 2005b) 



v(r) v 2 c (r) dr , 



(16) 



where u c (r) = ~\/GM(r)/r is the circular velocity profile, r a is the 
anisotropy radius, while the dimensionless kernel K is 



K(u)= -i / 1 



1 

^2 ' 



(17) 



for isotropic orbits (Tremaine et al. 1994; Prugniel & Simien 
1997) and 



f 1/2 / r / u a \ _, 



K(u, u a ) 



-sgn(u a -l)u a 



ui-l/2 



xC 



4 I U a U + 1 
U + U c 



("2-1/ 



(18) 



1 + -) cosh l u - - I -+' 
u 6\u 



M-l 
M+l 



(U a =l), 



where 



C-\x) 



cosh 



for u a > 1 



(Mamon & Lokas 2005b) for the anisotropy profile 



2 r + r a 



(19) 



(20) 



which Mamon & Lokas (2005b) (hereafter, ML) found to be a 
good fit to the anisotropy profiles of ACDM halos. 

Since the space density model enters equation (16) express- 
ing cr\ os (R) (through the tracer density v and the total mass M, 
which are related since we are considering single component 
mass models), we first need to determine the best fitting model 



to the distribution of particle radii in the stacked halo: we per- 
formed MLE of the NFW and Einasto models to the distribution 
of 3D radii of our stacked cluster. The minimum radius was cho- 
sen as 0.03 r2oo to avoid smaller radii, since our halo centers ap- 
pear to be uncertain to about 1% of the virial radius. We varied 
the outer radii of the fit, starting at r2oo- Since we will later fit 
the surface density profile out to the virial radius and beyond, 
we need to remember that the space radii extend beyond the 
maximum projected radius of the future surface density fits. So 
we also performed 3D fits beyond the virial radius: at 1.35 r2oo 
(which corresponds to the radius where A = 100, i.e. the largest 
radius where the halos should be close to virial equilibrium), and 
3 r2oo for a broader view of halos far beyond r 2 oo- 

When Prada et al. (2006) fit the density profiles of ACDM 
halos out to 2.7 r2oo, they found them to be well approximated 
by the sum of an Einasto model and a constant term pb? = QmPc- 
We therefore also experimented with the addition of a constant 
background component of density equal to the density of the 
Universe. In virial units, this background is expected to be equal 
to v bg = 3 Q m /(4;rA) = 3.6 x 10~ 4 (with Q m = 0.3 and A = 200). 



Table 1. MLE fits to the radial distribution of the stacked halo 



Model 



100 v b . 



NFW 


1 





4.08±0.05±0.17 


0.0021 


NFW 


1 


(0.036) 


4.10±0.05±0.17 


0.0023 


NFW 


1 


0.416 


4.26±0.07±0.23 


0.0069 


Einasto 


1 





4.00±0.05±0.17 


0.31 


Einasto 


1 


(0.036) 


4.02±0.05±0.17 


0.3 


Einasto 


1 


0.002 


4.00±0.07±0.20 


0.31 


NFW 


1.35 





4.20±0.05±0.14 


0.0023 


NFW 


1.35 


(0.036) 


4.23±0.05±0.14 


0.0021 


NFW 


1.35 


0.002 


4.20±0.04±0.15 


0.002 


Einasto 


1.35 





4.14±0.04±0.13 


0.016 


Einasto 


1.35 


(0.036) 


4.17±0.04±0.13 


0.0058 


Einasto 


1.35 


0.002 


4.14±0.04±0.14 


0.015 


NFW 


3 





4.80±0.04±0.24 





NFW 


3 


(0.036) 


5.01±0.04+0.24 





NFW 


3 


0.001 


4.80±0.03±0.19 





Einasto 


3 





4.26±0.03±0.22 


<10 15 


Einasto 


3 


(0.036) 


4.50±0.03±0.23 


<io- 10 


Einasto 


3 


0.042 


4.54±0.04±0.13 


<10" 



Notes. The Einasto models are for index m = 5. Column (2) is the 
maximum radius for the fits (the minimum radius is set to 0.03 ?'2oo)- 
Columns (3) and (4) are respectively the background density (multi- 
plied by 100, fixed if in parentheses) and concentration, while column 
(5) is the Kolmogorov-Smirnov test probability that the model is consis- 
tent with the distribution of radii. Virial units (with A = 200) are used 
for r lmx and v bg (r 2m and AW4)0' respectively), while c = r 20 o/r-2, 
where r_ 2 is the radius of density slope equal to -2. The outer radius 
of 1.35 r 2 oo corresponds to the virial radius, r l00 . Errors are 1 a: the first 
are the statistical errors from likelihood ratios and the second are from 
cosmic variance using 100 cluster bootstraps. Concentrations in bold 
face are those providing the highest P KS tests for given r max . 



Table 1 shows the resulting best-fit concentrations 9 obtained 
by MLE fits of NFW and m=5 Einasto models, plus an optional 
fixed or free constant background, to the distribution of 3D radii 
of the stacked halo. 10 The best-fit concentrations increase with 
the maximum allowed projected radius when the NFW model is 
used. This indicates the inadequacy of the NFW model at large 



9 In this paper, concentrations refer to rzoo/r-2- 

10 We also experimented with free index Einasto models: we generally 
found that the best fit index was in the range 4.6 < m < 5.2. 
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radii, as it fails to capture the steepening of the slope of the den- 
sity profile beyond the virial radius (Navarro et al. 2004). In fact, 
a Kolmogorov-Smirnov test (last column in Table 1) indicates 
that the NFW model is not an adequate representation of the dis- 
tribution of radii, whether a constant background is added or not, 
regardless of the maximum radius used in the fit. 

On the other hand, the concentration of the Einasto model 
appears to be somewhat less dependent of the outer radius, as 
also noted by Gao et al. (2008). At 0.03 < r/r 2Q0 < 1, the Einasto 
model is an adequate representation of the distribution of radii. 
At 0.03 < r < rioo = 1.35, the m-5 Einasto model is incon- 
sistent with the distribution of radii, but not by a large amount. 
However, the distribution of radii extending far beyond the virial 
radius, 0.03 < r/r2oo < 3 is not consistent with either NFW or 
m-5 Einasto models. 11 

The inclusion of the background in the fits leads to even 
higher concentrations when the maximum radius considered is 
3 r 200 . 



0.05 0.1 




1 



o 

3 0.9 
cc 

^ 0.8 h 
0.7 
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_L 
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+ Einasto 1 no 
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Fig. 10. Top: Space density profile (multiplied by r 2 ) of the stack of the 
93 halos, with best maximum likelihood fits for r/r 2 oo from 0.03 to 1 
(solid), and 3 (dotted without background, dashed curves with best-fit 
background, see Table 1) for the NFW (red) and Einasto (blue) mod- 
els. The best fit NFW models to 3 r 2 oo with and without background are 
indistinguishable (see Table 1). The errors are from 100 cluster boot- 
straps. Bottom: Ratios of measured to fit densities. 



11 For r max = 3 r 2 oo, the KS test showed that the free m Einasto model 
with free or fixed background (or without any) failed to provide an ad- 
equate representation of the distribution of radii. 



The density profile of the 93 stacked clusters is shown in 
Figure 10 for maximum fit radii of r max = 1 and 3 r^o- The max- 
imum likelihood NFW model produces clearly worse fits to the 
density profile than the maximum likelihood Einasto model for 
r < 3 r 2 oo, while the NFW model reproduces better the mea- 
sured density profile at r — 3 r 2 QQ. As clearly seen in the bottom 
panel of Figure 10, neither model is adequate near 2 r 2 QQ, even 
when considering the cosmic variance measured by our cluster 
bootstraps (see the error bars in the top panel of Fig, 10). 

In most of what follows, we restrict our analysis to R < r 2 oo. 
For these analyses, we adopt the c — 4, NFW and m-5 Einasto 
models, as these models are simple and the latter is consistent 
with the radial distribution without and with a constant back- 
ground. 



1.4 



1.2 




lL 



0.4 



0.2 p 



-0.2 



0.01 



0.1 

r / r 200 



Fig. 11. Velocity anisotropy profile (including streaming motions: 
eqs. [6] and [7]) of the stack of the 93 halos. The error bars are from 
100 bootstraps on the 93 halos. The curve shows the weighted^ 2 fit (in 
the range 0.03 < r/r 2 oo < 1) of the Mamon-Lokas anisotropy (2005b) 
model (eq. [20]) with r a = 0.27 r 2 oQ (the solid portion of the curve high- 
lights the region where the fit was performed). The purple dashed hori- 
zontal line indicates the fully isotropic case. 



The radial profile of velocity anisotropy (eq. [6]) of the 
stacked cluster, shown in Figure 11, is reasonably well fit by 
(reduced^ 2 = 1.4) the ML model 12 (eq. [20]), with anisotropy 
radius (where /3 reaches its half-value of 1/4) r a = 0.275 ± 
0.020r2oo, as found with a weighted x 1 fit for logr/r2oo be- 
tween —1.5 and 0. With a concentration parameter c - 4.0, 
the anisotropy radius is r a 1.1 r_2, i-e. the velocity anisotropy 
reaches its value intermediate between the center and the outer 
regions close to the radius of density slope —2. 

We now measure the velocity dispersion of the stacked pro- 
file using different schemes for interloper removal to find a 
scheme that produces a velocity dispersion profile (on the data 
with Hubble flow and the velocity cut) in agreement with the pre- 



12 Other anisotropy models such as constant and Osipkov-Merritt 
(Osipkov 1979; Merritt 1985) produce much worse best fits (reduced 
X 1 — 5 and 20, respectively). 
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dictions (with no Hubble flow nor velocity cut) for the c — 4.0 
NFW and Einasto models whose density profiles we just fit in 
3D. Since the surface density profile enters equation (16), and 
no analytical formula is known for the Einasto model, we de- 
rived an accurate approximation for the Einasto surface density 
profile (for a large range of projected radii and of indices rri) in 
Appendix A (eqs. [A. 15] and [A. 14]). 

The red open triangles in the top panel of Figure 12 show the 
line-of-sight velocity dispersion profiles "measured" with our 
standard velocity cut at 3 times the predicted isotropic line-of- 
sight velocity dispersion (hereafter <x™) for an NFW model with 
concentration c - 4.0 (as measured in 3D, see above). In com- 
parison, cr™ (red solid curve in the top panel of Fig. 12) is typi- 
cally 10% lower than the "measured" velocity dispersion profile 
for radii R < 0. 1 7-200- This discrepancy is decreased to 4% if one 
compares the measured velocity dispersions after clipping at 3 
times the line-of-sight velocity dispersion, computed with the 
ML anisotropy (hereafter crj^, eqs. [16] and [18], black filled 
triangles) to crj^ (black solid curve in the top panel of Fig. 12). 
This suggests that the 3 <x clipping generally used is too liberal. 
A near perfect match (typically better than 1 % for R < 0.8 r2oo) 
is obtained by cutting at 2.7 crM?* (black filled circles vs. black 
solid curve in the top panel of Fig. 12). 

When the m-5 Einasto model is used to compute cr\ os (R) be- 
fore applying the velocity cut, the best match between the mea- 
sured and predicted line-of-sight velocity dispersion profiles is 
for a cut at k — 2.6 (bottom panel of Fig. 12). 

The Hubble flow (HF) causes a shallower slope at projected 
radii close to the virial radius (one notices in both panels of 
Fig. 12 an inflection point in the measured profiles [filled circles] 
of log <xi os vs. log R near half a virial radius). Indeed, we obtained 
results similar to those of Figure 12 when we did not incorporate 
the HF to the peculiar velocities of the simulation: the measured 
c\os(R) fell more sharply, with no inflection point, even some- 
what more sharply than predicted by the Einasto model (because 
the velocity anisotropy without the HF is more radial at 4 r2oo 
in comparison with the case where the HF is incorporated, where 
4 f"2oo roughly corresponds to the turnaround radius where the ve- 
locities are mostly tangential). So, although the steeper Einasto 
density profile ought to catch better the steeper line-of-sight ve- 
locity dispersion profile at large projected radii, the NFW model 
performs slightly better, because its shallower line-of-sight pro- 
file mimics better the effects of the Hubble flow. 

In summary, Figure 12 indicates that if one wishes to recover 
the correct line-of-sight velocity dispersion profile, one should 
use 2.6 or 2.7 cr clipping instead of 3 o~ clipping, where the line- 
of-sight velocity dispersion is either measured or modeled with 
anisotropic velocities. 

The choice of model and k is not obvious. We prefer the 
NFW model, as it is simpler and, with k = 2.7, it presents 
a slightly better match between measured and predicted line- 
of-sight velocity dispersion profiles than does the m-5 Einasto 
model (compare the ratios of measured to predicted o~\ os (R) in 
both plots of Fig. 12, especially at large radii). 

Mass modelers of clusters may wish to avoid performing the 
integral of equation (16) with the kernel of equations (18) and 
(19). The line-of-sight velocity dispersion profile (eq. [16]) for 
the NFW and m-5 Einasto models with ML anisotropy with r a - 
r_2 can be approximated as 
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Fig. 12. Line-of-sight velocity dispersion profiles of the stacked virial 
cone, cutting at 3 (triangles) or 2.7 (circles) o~ los (K), assuming the best 
fitting (c = 4.0) NFW (top) or m = 5 Einasto (bottom) model with 
isotropic velocities (red open symbols) or slightly radial ML anisotropy 
(eq. [20]) with r a = r_2 (black filled symbols). The error bars are from 
100 bootstraps on the particles within each bin of projected radii. For 
clarity, the isotropic and ML symbols are shifted by 0.01 dex left- 
wards and rightwards, respectively. The curves show the predicted line- 
of-sight velocity dispersions (for no Hubble flow, eq. [16]) assuming 
isotropy (red curve, eq. [17]) or ML anisotropy with r a = r_ 2 (black 
curve, using eq. [18]). The ratios of measured to predicted cr los (i?) are 
shown in the lower frames of each plot (same colors and symbols as 
upper frames). 



0-l 0s (ft) 



VGM(r_ 2 )/r_ 



dex 



logi 



R 

r-2 



where the coefficients are given in Table 2 for both models. 
(21) These two approximations are accurate to better than 0.5 % (rms) 
for 0.0032 < R/r_ 2 < 32. Then, one can write cr los (/?) in terms 
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Table 2. Coefficients for cn os (7?) approximation (eq. [21]) 





NFW 


m=5 Einasto 




— U. 14 / o 


— U. 1 JZU 


"I 


-0.1109 


-0.1242 


a 2 


-0.1357 


-0.1637 


"3 


0.001948 


-0.01688 


"4 


0.02317 


0.01892 


"5 


0.0006310 


0.001844 


06 


-0.003234 


-0.002044 


a-j 


-0.0006370 


-0.0004103 



of v v using equation (21) and 

ln(c + 1) - c/(c + 1) 
(ln2- 1/2) c 

y(3m,2mc 1/m ) 



GM(r- 2 )/r-2 



(NFW) , 



(22) 



y(3m, 2m) c 



(Einasto) 



(e.g., Navarro et al. 1996 for NFW and trivially derived from the 
Einasto mass profile first derived by Mamon & Lokas 2005a). 

In the absence of information on the mass profile (e.g. from 
X ray observations), neither the concentration parameter, c, nor 
the scale radius r_2 are known, so one has to work iteratively, 
first guessing a plausible value of c, applying the velocity filter, 
then re-estimating c from the data and re-applying the velocity 
filter. This process should converge in a one or two iterations. 

6. Interloper statistics after the velocity cut 

We now show the statistics of interlopers after our adopted ve- 
locity cut. The motivation is to allow observers to compare with 
their own data. With the velocity cut (k = 2.7), the line-of-sight 
distances are now effectively limited to +17 r^oo from the cen- 
ter of the stacked halo (eq. [5]). The qualitative features of the 
projected phase space distribution are robust to variations of the 
method to cut the velocities. 

Figure 3 shows the velocity cut at +2.7 cri os (R) — with our 
adopted NFW model with ML anisotropy with anisotropy radius 
r a = r_2 — on top of the projected phase space (using eqs. [21] 
and [22]). Only 0.4% of the halo particles are rejected by the 
2.7 cr velocity cut, which is a low enough fraction that the shot 
noise in the structural and kinematical modeling is not signifi- 
cantly increased. The velocity cut in Figure 3 seems very rea- 
sonable as it is close to optimizing the completeness of the se- 
lection of particles within the virial sphere. However, less than 
17% of the interlopers are identified as such by the velocity cut. 
Therefore, the great majority of interlopers cannot be removed 
by a velocity cut. 

The green curves in Fig. 4 show the velocity cut on top of 
the phase space density distribution. The velocity cut appears 
to occur in a region where the interloper phase space density is 
roughly constant. 

This can be seen in a clearer fashion in Figure 6, where the 
velocity cut is shown as green vertical lines. While the high- 
est velocity interlopers are removed by the velocity cut, there 
remains signs of the field component, which we had identified 
with particles beyond 8 virial radii, at low (R < 0.4 r2oo) pro- 
jected radii. 

The fraction of particles outside the virial sphere is displayed 
in Figure 13. Interestingly, at R > 0.8 rioo (magenta histogram) 
the interlopers account for over 60% of all particles, regardless 
of the particle velocity up to the velocity cut (filled circles). But 



m 
u 

<D 

a. 

o 




Fig. 13. Fraction of interlopers as a function of line-of-sight velocity for 
all projected radii R < f2oo (thick black histogram) and in bins of pro- 
jected radius: R/r 2Q0 = 0-0.2, 0.2-0.4, 0.4-0.6, 0.6-0.8, and 0.8-1 (thin 
histograms), increasing upwards. The filled circles show the 2.7 a\ 0S (R) 
(from eqs. [21] and [22]) velocity cut for the c=4 NFW model with 
r a =r_2 ML anisotropy (eq. [20]). 



even at smaller radii, 0.4 < R/r2oo < 0.6, interlopers account for 
over 20% of all particles again for all velocities up to the cut. So, 
unless one limits one's kinematical analysis to very small cluster 
apertures, one cannot avoid being significantly contaminated by 
interlopers. 

While there is no gap in the velocity distribution of parti- 
cles (Fig. 6), Figure 13 shows local minima of the interloper 
fraction for all bins of projected radii, except the outermost one. 
Regardless of the application of a velocity cut, these local min- 
ima occur at lower velocity (1.3 to 1.5 u v ) than the inflection 
points of the interloper density in projected phase space (1.6 to 
2.6 1\, as seen in Fig. 6). The local minima occur at velocities that 
decrease with projected radius, suggesting that our local kcti os 
cut is preferable to a global one, since cr\ os (R) decreases with 
R for R > 0.1 r2oo (see Fig. 12). These local minima arise be- 
cause the interloper system has a lower velocity dispersion than 
the halo system: o~\ = 0.58 (eq. [10]) while after the velocity cut 
the aperture velocity dispersion of the global stacked virial cone 
(thus including both halo particles and interlopers) is 77 = 0.65, 
i.e. 5% higher than predicted by Mauduit & Mamon (2007) for 
an isotropic NFW model, which is not surprising given the radial 
anisotropy of the halos (Fig. 11). 

The surface density profile of the stacked halo is shown in 
Figure 14. The surface density profile of the interlopers is flat 
with small fluctuations around the mean values = 0.114 and 
0.096 N v r~ 2 , measured in the stacked virial cone, respectively 
before and after the velocity cut. 13 In comparison, our model of 
the surface density of interlopers (eq. [12]) combined with our 



13 Figure 14 shows interloper surface densities that are lower, at R < 
0.7 r 2 m, man 0.114 and 0.096 ./V v r~ 2 , respectively before and after the 
velocity cut, but most of the particles lie within the highest bins of log 
projected radius. 
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MLE values for A, B and cr ; yields mean interloper densities of 
0. 1 14 and 0.094 N v r~ 2 , respectively before and after the k = 2.7 
(k — 1 .76 with 77 = 0.65) velocity cut. The general agreement is 
excellent. 
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Fig. 14. Top panel: Surface density profile of global stacked cone (black 
histogram, raised up by 0.02 dex for clarity), as well as halo members 
(r < r 2m , red histogram) and interlopers (r > r 2 oo) before (thin) and 
after (thick blue histograms) the velocity cut. Poisson errors are only 
shown for the interlopers after the velocity cut for R < r 200 . Also shown 
(curves) are maximum likelihood m=5 Einasto model fits (after the ve- 
locity cut), in the range 0.03 r 2 oo to 1 (purple) or 3 (brown and green) 
r 2 oo, without (purple and brown) or with (green) an additional free con- 
stant background component (dashed green line). Bottom panel: Ratios 
of measured to fit surface densities (after the velocity cut). For clarity, 
Poisson errors are only shown if larger than the symbol size. 



Note that, at projected radii beyond the virial radius, all par- 
ticles are interlopers, so the surface density of interlopers is not 
constant but decreases, to first order, as the NFW or Einasto 
models. While the total surface density profiles of the popular 
NFW and Einasto models for ACDM halos are convex in log 
surface density vs. log projected radius, an important additional 
background term in the surface density would lead to an inflec- 
tion point and subsequent concavity at some radius. Such a fea- 
ture would lead to poor fits of single NFW or Einasto profiles. 
Now, within the virial radius, no such inflection point and outer 
concavity are seen in Fig. 14 for the total surface density pro- 
file. However, extending the surface density profile out to three 
virial radii, as illustrated in Figure 14, one does see the inflec- 
tion point of the total surface density profile (near 1 .6 7-200 after 
the velocity cut and right at 7-200 before). This points to an ad- 



ditional background of surface density. This requirement for the 
additional background component is confirmed by the fairly flat 
ratios of data over model (bottom panel of Fig. 14) for the case 
where the background is fit, in comparison with larger residuals 
for R > 7-200 for the fits without a background. 

Such a background is expected, since the density profiles of 
ACDM halos is the sum of an Einasto model and a constant 
background corresponding to the mean density of the Universe 
(Prada et al. 2006). Integrating along the line-of-sight (eq. [1]) 
within the sphere of radius r max (eq. [5]), one deduces that the 
surface density profile of (foreground/)background structures is 



2 bg (#) = 2 Q m p c ^jr^-R 2 . 



(23) 



Since r max » r v (eq. [5]), £bg is roughly constant for R 5 r v : 
^bg — 2 fi m p c y max , (24) 
which in dimensionless virial units (N v r^ 1 ) becomes 
a £b g 3 ?]Kn m 

^bg = ,, , n ;=- , (25) 



N v /r; n V8A 

using equation (3). With the velocity cut, k = 2.7 and equa- 
tion (25) yields t hg 0.0126 for fi m = 0.3, A = 200, and 
77 = 0.65 (see above). Without the velocity cut k = t]k = 4 and 
one obtains Zb g - 0.0286. According to equation (23), the rela- 
tive drops of 2b g from R = to R mdx = 1 or 3 r v are 0.2% and 
1.6%, respectively, so the approximation of a constant surface 
density background is adequate. 

This background corresponds to the velocity-independent 
component of the interloper surface phase space density, which 
in our model (eq. [8]) is the B term, which produces a mean sur- 
face density of t hg = 77 x 2.7 B = 0.0132 (with B = 0.0075 from 
eq. [10] and again 77 = 0.65). This agrees with the previous value 
to within 4%. This means that the constant field term in the ve- 
locity distribution corresponds precisely to the additional halos 
outside the test halo. In any event, the total surface density of a 
cosmological structure is the sum of the surface density of that 
structure and a constant background. 

Although the mean surface density of interlopers is roughly 
independent of halo mass (Fig. 9), it may vary from cluster 
to cluster. Figure 15 shows the surface density of each of the 
93 halos. As seen in Table 3, the arithmetic mean value of Ej 



Table 3. Statistics of interloper surface densities E ; (in N v r v 2 ) 



Velocity cut 


No 


Yes 


Arithmetic mean 


0.116 


0.098 


Geometric mean 


0.102 


0.085 


Standard deviation of log £j 


0.216 


0.229 


Spearman rank correlation 


0.025 


-0.008 


Probability 


0.34 


0.44 



Notes. Spearman rank correlation is between mean surface density of 
interlopers and halo mass. Probability is of having a stronger correlation 
by chance. 



matches well the value of the stacked virial cone, regardless of 
the velocity cut, which has only a minor effect on the statistics 
of interlopers (while the geometric means are 15% lower). 
But the dispersion in logE; is as high as 0.22, close to cr(logA) 
(eq. [14]), so that the relative dispersion of 2; is as high as a 
factor lO 22 * 5/3. 
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Fig. 15. Mean interloper surface density (in virial units) versus halo 
mass (for A = 200). The open black and filled blue circles show the 
93 halos each measured along 3 orthogonal viewing directions, before 
and after the velocity cut, respectively. The right panel provides the fre- 
quency of the mean interloper surface density (with logarithmic bins), 
before (dashed black) and after (solid blue histograms) the velocity cut. 



The fraction of interlopers in the stacked virial cone is 



fi = 



Ni 



71% 



N h +Ni l+7r£i 



(26) 



(where indices 'h' and 'F correspond to halo and interloper par- 
ticles, respectively), where the second equality of equation (26) 
made use of £h = %/[N v /rl] = l/n, by definition. This yields 
/i = 27% before and 23.1+0.1% after the velocity cut, where the 
error is both from binomial statistics and from a bootstrap on all 
the particles; a bootstrap on the halos leads to an error of 0.6%; 
propagating (with eq. [26]) the error on the mean of % (from the 
standard deviation of % given in Table 3) also leads to an error 
on fi of 0.6%; finally, the standard deviation of the interloper 
fraction for the three stacked cartesian virial cones is 1.7%. We 
adopt the error estimate of 0.6% on f\ for the later discussion. 

The small decrease in interloper fraction from before to af- 
ter the velocity cut confirms our finding that the large majority 
of interlopers have too low velocities to be filtered by velocity. 
Cosmic variance causes huge fluctuations in the fraction or sur- 
face density of interlopers (2/3 of the mean value, independent 
of the presence of a velocity cut), with roughly a log-normal dis- 
tribution (see the right panel of Fig. 15 and Table 3). The last two 
lines of Table 3 indicate that there is no statistically significant 
correlation of surface density of interlopers with halo mass. 



7. Biases in concentration and anisotropy? 

What are the effects of the Hubble flow on estimates that ob- 
servers make on halos, e.g. the concentration and the velocity 
anisotropy of the distribution of their tracer constituents (i.e. 
galaxies in clusters)? 



7.1. Effects of the Hubble flow 

We begin by a naive comparison of the observable distributions 
with and without the Hubble flow, before comparing the concen- 
tration and anisotropy measured by an observer with the corre- 
sponding quantities we directly infer in 3D from the cosmolog- 
ical simulations. Admittedly, the distribution of velocities with- 
out the Hubble flow is not fully realistic, since the cosmological 
simulation solved equations for comoving coordinates in an ex- 
panding universe. 14 
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Fig. 16. Difference of velocity moments with Hubble flow and veloc- 
ity cut (HF) and without Hubble flow or velocity cut (noHF): log sur- 
face density (dotted black), log line-of-sight velocity dispersion (solid 
red), and line-of-sight velocity kurtosis (dashed blue). The error bars 
are based upon Poisson errors for the surface density and bootstraps 
within the radial bin for the log dispersion and the kurtosis, where the 
error on the difference is the square root of the sum of the square errors. 



Figure 16 shows the changes in the radial profiles of surface 
density and line-of-sight velocity dispersion and kurtosis, 15 once 
the Hubble flow is added to the peculiar velocities. One striking 
feature of Fig. 16 is that the Hubble flow leads to a lower sur- 
face density profile at large radii. This cannot be a consequence 
of the restriction of the line-of-sight of the halo component to 
±19r2oo (see Sect. 1), because the NFW surface density with 
line-of-sight limited to the sphere of that radius (appendix B) 
matches the NFW surface density projected to infinity to better 
than 1.4% relative accuracy for R < 3 r2oo (for c = 4). Instead, 
it is the integral along the line of sight of the constant den- 
sity (foreground/background) component that diverges when no 
Hubble flow is present, and is limited to half a box size here: 



14 In fact, static universes are never simulated in a cosmological 
context, because of their lack of realism, given the expansion of the 
Universe as seen in the Hubble law, and also because of the lack of 
knowledge of suitable initial conditions in such a static universe. 

15 The reader should not confuse the line-of-sight velocity kurtosis k ]os 
with the velocity cutoff in units of line-of-sight velocity dispersion (k) 
or of virial velocity (£). 
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L/2 - 96 h Mpc, which corresponds to roughly 100 virial 
radii. In any event, the lower (40% lower at r2oo) surface den- 
sity profile found when the Hubble flow is added to the peculiar 
velocities might explain the lack of concavity in the (log-log) 
surface density profile within the virial radius (Fig. 14). 

We noticed in Sect. 5 that the line-of-sight velocity disper- 
sion profile showed an excess at radii near the virial radius. 
Figure 16 confirms that the line-of-sight velocity dispersion pro- 
file is gradually overestimated at large radii, while the surface 
density profile is much more biased beyond half a virial radius, 
being underestimated at large radii. The effects of the Hubble 
flow on the surface density and line-of-sight velocity dispersion 
profile are both small at very low projected radii. 

Finally, our sharp cut (Sect. 5) in the distribution of line-of- 
sight velocities when the Hubble flow is incorporated implies 
that the line-of-sight velocity kurtosis is underestimated ( espe- 
cially at large projected radii). 

7.2. Concentration 

Does the excess of interlopers at large projected radii lead to 
lower values of the concentration parameter in the fits of the pro- 
jected NFW and Einasto profiles to the surface density profiles 
of clusters? 

Table 4 shows the MLE fits (see appendix C) of the NFW 
and m-5 Einasto surface density profiles to the distribution of 
projected radii of the three cartesian stacked cones. Different fits 
were performed with variations in the maximum allowed pro- 
jected radius, R max , the presence of a constant (fixed or free) 
background term, and the possible removal of high-velocity out- 
liers. The NFW surface density and projected number (or equiv- 
alently projected mass) profiles, required for the normalization 
of the probability used in the MLE, are given by Bartelmann 
(1996) and, in another form by Lokas & Mamon (2001). The sur- 
face density and projected number (mass) profiles of the Einasto 
model are not known in analytical form, so we have derived ac- 
curate approximations in appendix A. For these MLE, we adopt 
equation (A. 17) with equation (A. 6) for the surface density pro- 
file and equation (A. 8) with equation (A. 6) for the projected 
mass profile. When a constant background term is included in 
the fits, it is either free (column 8) or fixed at £bg —kB — 0.0286 
and 0.0126 without (k = 4) and with (k — tjk — 1.75) the veloc- 
ity cut (see Sect. 6). 

The concentrations measured on the projected radii with sin- 
gle component fits are always smaller than the characteristic 
value found in 3D (c = 4.0, 4.1, and 4.5, see bold values in 
Table 1). The best-fit concentrations are very low when the fits 
are performed out to 3 r2oo (unless a velocity cut is performed 
or background component added to the model). Note that when 
the background is fitted together with the concentration and no 
velocity cut is performed, the best-fit value for the background 
can be over a factor two off from the value expected from equa- 
tion (24), even with maximum projected radii of 3 r2oo- The KS 
tests indicate that for most combinations of maximum projected 
radius, presence or absence of the velocity cut and how the back- 
ground is handled, the m-5 Einasto model usually provides a 
better representation of the distribution of projected radii than 
does the NFW model. Finally, the errors in Table 4 indicate that 
the cosmic variance of stacks of 93 clusters, measured using the 
standard deviation of the three cartesian stacked cones with the 
gapper 16 estimate of dispersion, which is most robust to small 

16 The gapper dispersion of a vector x of length n is s = ^/K/[n(n - 
1)] 2Xi' i(n - i) (Xj+i - x t ) (Wainer & Thissen 1976). 



sample sizes (Beers et al. 1990), are much greater than the in- 
trinsic fitting errors. 

Table 5 shows the bias in concentrations, where the bias is 
the ratio between the concentration measured with the 2D fit 
(Table 4) and the best-fitting of the concentrations found with 
the 3D fits (highlighted in bold in Table 1), for given r max /r2oo 
and different models and backgrounds. The concentration pa- 
rameters found in the fits of the surface density profile are un- 
derestimated by typically 16%+7% when we make no velocity 
cut, limit the projected radii to t-2oo> an d do not incorporate a 
constant background in the fit. This underestimate of the con- 
centration is statistically significant (marginally so for NFW). 
The concentration bias gets worse as we increase the maximum 
projected radius: the concentration is underestimated by 1/4 at 
R < rioo and by as much as 2/3 with R < 3 ^200- 

But when we make the velocity cut at 2.7 (NFW) or 2.6 
(Einasto) cr\ os (R), where cr\ os (R) is obtained from equations (21) 
and (22), the biases in the concentration parameters are typically 
reduced by a factor two. When we limit the analysis at R < r2oo, 
the concentration parameters found in the fits of the surface den- 
sity profiles are (only) underestimated by typically 6+6%. This 
low bias is no longer statistically significant given our fit and 
cosmic variance errors. 17 However, extending the analysis to 
R = fioo, the concentration is still biased low by as much as 
11% (which is marginally significant), despite the velocity cut. 
And if we go all the way to 3 ^200, the bias is still very strong, as 
the concentration is underestimated by over 1/3. 

One may wonder whether one can recover the concentration 
parameter more accurately with two-component fits to the set 
of projected radii than with the single-component fit, since we 
found an additional background needs to be added (Sect. 6). For 
example, Lin et al. (2004) fit a projected NFW model plus a 
constant surface density background (hereafter NFW+bg), with 
no velocity cut, for projected radii 0.02 < /^max/^oo < 2.5. We 
tested the single and double-component models in the optimistic 
case of our stacked cluster with nearly 3 x 10 5 particles. As seen 
in Table 5, the 3D concentration parameter is recovered better 
by the two-component models, regardless of the velocity cut, al- 
though the improvement is not always statistically significant. 
Note that the background is well recovered with R mdx = 3 r2oo 
and a velocity cut for both the NFW+bg and Einasto+bg models 
(Table 4). 

In summary, the concentrations measured in 2D recover best 
the 3D values once the velocities are filtered and especially once 
a background is included in the fit ( even when limited to the 
virial radius). 

Note also that one should not attempt to model the surface 
density profile as the sum of the halo term (with line-of-sight 
limited to the sphere, with the formulae of appendix B for the 
NFW model) and a constant background, because the total sur- 
face density profile decreases smoothly beyond the virial radius 
in ways that are not simple to model, for example with a spher- 
ical halo and a constant background. Moreover, for fits where 
the projected radii are limited to the virial radius, these spherical 
plus background fits are not recommended because the back- 
ground is not constant but rises with radius (Figs. 5 and 8) and 
is less easy to model than the surface density with line-of-sight 
integrated to infinity. 



17 We also experimented with concentration fits on the projected radii 
of all particles within 3 cr\ os (R) instead of 2.7, but the changes were very 
small (less than 1%, with slightly worse underestimates of the concen- 
tration for the single component fits). 
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Table 4. MLE fits to the distribution of projected radii of the three cartesian stacked cones 
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Notes. Column 1: model (NFW or m=5 Einasto); column 2 (i? max ): projected radius of the cone in which the stacked cluster is built, in units of 
r 2 go', column 3 (u-cut): presence (Y) or absence of the velocity cut with k = 2.7 (NFW) or 2.6 (Einasto); columns 4-5, 6-7, 9-10: mean best-fit 
concentration (c = r200/V-2) from projected radii and probability that distribution of projected radii is consistent with model using a Kolmogorov- 
Smirnov test (PrsX f° r fits without a background (cols. 4 and 5), with a fixed background (L bs , = 0.0286 [no velocity cut] or 0.0126 [with velocity 
cut], cols. 6 and 7) or a free background (cols. 9 and 10), with best-fit value (100 times the geometric mean and error on its logarithm in parentheses) 
given in column 8. The minimum projected radius is set to 0.03 r 2 oo- For the Einasto model, we adopt the approximation to the surface density and 
projected number (mass) profiles given in appendix A. The errors on c are statistical (first) and a measure of the cosmic variance term estimated 
by the gapper (Wainer & Thissen 1976, see Beers et al. 1990) standard deviation of the MLE values for the three projection axes. 



Table 5. Concentration bias of 2D fits 
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NFW 




1.35 
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Notes. The biases highlighted in bold (respectively blue italics) show the cases where the best-fit surface density profile was consistent with the 
data to better than 5% (between 0.01% and 5%) confidence (see last column of Table 4). The errors are the statistical (MLE fit) and cosmic variance 
(from the 3 cartesian stacked halos) errors of Table 4 added in quadrature together and with the analogous errors on the best-fit (bold in Table 1) 

3D concentration, with or 2 (c 1D /c b ^) = cr 2 (c 2D )/ (c^f + o- 2 (c^ st ) <c 2D > 2 / (c^f . 



7.3. Velocity anisotropy 

The signature of the Hubble flow on the shape of the line-of-sight 
velocity dispersion profile (Fig. 12) suggests that the velocity 
anisotropy that is recovered may be affected. Figure 17 shows 
the result of the non-parametric anisotropy inversion (first de- 
veloped by Binney & Mamon 1982, but we use here the simpler 
algorithm by Solanes & Salvador-Sole 1990), which computes 
the anisotropy profile assuming a smooth representation of the 
line-of-sight velocity dispersion profile and a mass model. Here, 
we adopt a c=4 NFW model and fit polynomials to the binned 
logcr los vs. log R. 18 

In the region where the order of the polynomial fit does not 
matter (0.1 < r/rioo < 1), the recovered anisotropy profile re- 
produces very well the one measured in three dimensions (points 
in Fig. 17), although, beyond 0.2r2oo, the recovered anisotropy 

18 We cannot employ the analytical approximation of equation (21) 
to the line-of-sight velocity dispersion profile for an NFW model with 
r a = r_2 ML anisotropy, because we place ourselves in the context of an 
observer who wishes to measure the velocity anisotropy with no prior 
on it: (s)he is thus forced to use a smooth representation of the observed 
line-of-sight velocity dispersion profile. 



profiles are slightly more radial than that measured in 3D. This 
bias towards more radial motions appears statistically signifi- 
cant, since in all six radial bins where there is a ^ 1 <x offset, 
this offset is in the same direction (probability of 2~ 5 = 3%). 
Therefore, the Hubble flow produces only a slight radial veloc- 
ity anisotropy bias in the envelopes of halos. 

8. Summary & discussion 

This work analyzes the distribution of particles in projected 
phase space around dark matter halos in cosmological simula- 
tions. The particles are split among halo particles within the 
virial sphere and interlopers within the virial cone but outside 
the virial sphere (Fig. 2). The reader should be careful that the 
analyses presented here cannot be directly applied to observa- 
tions of clusters of galaxies, as they work with halo particles in- 
stead of galaxies within clusters, and assume the halo centers to 
be determined quite precisely (from real space measurements). 

We find a universal distribution of interlopers in projected 
phase space, i.e. with little dependence on halo mass (Figs. 7c 
and 9). In particular, we note that velocity cuts cannot distinguish 
the quarter of particles that are interlopers from those in the virial 
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0.1 1 

r / r 200 

Fig. 17. Velocity anisotropy profiles (including streaming motions: 
eqs. [6] and [7]) of the stack of the 93 halos. The points are the mea- 
sured velocity anisotropy (same as in Fig. 1 1, again with uncertainties 
from 100 bootstraps on the 93 halos) and the curves are recovered from 
anisotropy inversion assuming the c=4 NFW model (solid curves), for 
three polynomial fits (orders 2, 3, and 4 in log-log space) to the mea- 
sured line-of-sight velocity dispersion profile (after the k = 2.7 ve- 
locity cut using the NFW model with ML anisotropy: solid dark red, 
dashed green, and dotted blue for orders 2 to 4, respectively). The pur- 
ple dashed horizontal line indicates the fully isotropic case. Note that 
the 4th order polynomial (blue) extrapolates poorly the line-of-sight ve- 
locity dispersion profile at very low and very high projected radii. 



sphere (Fig. 13), as was previously noted by Cen (1997). We find 
that the distribution of interlopers in projected phase space dis- 
plays a roughly constant surface density (Figs. 5 and 8) and a 
distribution of line-of-sight velocities that is the sum of a quasi- 
gaussian component, caused by the halo outskirts (out to typi- 
cally 8 virial radii, Fig. 7b) and a uniform component caused by 
particles at further distances from the halo (Figs. 6 and 7). 

The cosmological simulations allow us to optimize the ra- 
tio of maximum velocity to line-of-sight velocity dispersion that 
recovers the latter quantity. Although this may seem to be a cir- 
cular argument (since a-\ m (R) depends on the velocity cut), it has 
been widely used in the past, usually in iterative form, with a 
3 cr cutoff. We find that this cutoff is not restrictive enough and 
causes an overestimate of the line-of-sight velocity dispersion 
profile (based upon mass and velocity anisotropy models derived 
from the cosmological simulations): up to 10% for the isotropic 
NFW velocity cut, which is reduced to 5% for the ML anisotropy 
velocity cut (Fig. 12). We recommend instead a velocity cut at 
2.7 <Ti os (/?) on the best iterative fit to the line-of-sight velocity 
dispersion for the NFW model with r a = r_2 ML anisotropy. 
Alternatively, one can use a velocity cut at a- = 2.6 for the m=5 
Einasto model, modeled (again) with r a = r_2 ML anisotropy, 
but this underestimates the line-of-sight velocity dispersion near 
the virial radius (Fig. 12). 

We illustrate (Figs. 3, 4, 6, 13, and 14) how the distribution of 
particles in projected phase space is altered once the high veloc- 
ity interlopers are rejected with this new velocity filter (besides 



limiting the line-of-sight to typically + 17r2oo, the main effect 
is to remove the flat velocity component). The fraction of inter- 
lopers within the virial cone drops from 27% (with an observer 
at distance D = 90 /i" 1 Mpc) to 23.1 ±0.6% (independent of D 
for D ^ 17 (r2oo)) when the velocity cut is applied (where the 
uncertainty is taken from the end of Sect. 6). 

This fraction of interlopers can be directly inferred from the 
NFW or Einasto model 

= M p (r 20 o)-l+7ri: bg 
M p (r 2 oo) + 7rt bg 

where M p = Mp/Mooo is the projected virial mass in virial units 
(i.e. in units of the mass within the virial sphere), while Zbg is 
given in equation (25). For the NFW model, one then obtains fi 
= 26.8% and 24.0%, respectively before (E bg = 0.0286) and af- 
ter (£bg = 0.0126) the velocity cut, while with the m-5 Einasto 
model with c — 4, the corresponding percentages of interlopers 
are 25.7% and 22.8%. These theoretical predictions are in excel- 
lent agreement with the fractions obtained from the simulations. 
Note that the omission of the background (2bg) term in equa- 
tion (27) reduces fi by typically 10% in relative terms, relative 
to the fractions after the velocity cut. Applying equation (27) to 
models of different concentrations leads to roughly a power-law 
variation of f with slope -0.32 (NFW) or -0.49 (m-5 Einasto). 
Therefore, the fraction of interlopers should be (slightly) more 
important in the more massive halos, since they have (slightly) 
lower concentrations (Navarro et al. 1997; Maccio et al. 2008). 

In comparison, using 62 clusters from the same simulation 
as the one we have analyzed (we have 53 clusters in common), 
Biviano et al. (2006) found that among particles selected in 
cones of projected radius 1.5 /i -1 Mpc around cluster-mass ha- 
los (after their velocity cut), 18+1.4% of them lie outside the 
sphere of the same radius. 19 Their halos have a median virial ra- 
dius of f-200 = 0.93 hT l Mpc (1.08 times our median) and hence a 
virial mass of M200 = 1.9 x 10 l4 h~ l M G . Their Figure 7 indicates 
that their velocity cut is roughly 1180, 1105, and 780 kms -1 , 
at projected radii 0.6, 1.0 and 1.5 lr l Mpc, respectively. Since 
NFW concentration scales as M~ ()1 (Navarro et al. 1997; Maccio 
et al. 2008), their median concentration should be 3.9, hence 
their scale radius should be 930/3.9 = 238/i _1 kpc. Assuming 
an NFW model, we deduce that their median circular velocity 
at the scale radius is 914 km s _1 , and find that their velocity cuts 
correspond to k — 2.1,2.2, and 1.8, at the three projected radii 
chosen above. These fractions are consistent with the values of 
k one can read off of Figure 3 of Wojtak et al. (2007) that illus- 
trates the same velocity cut model (den Hartog & Katgert 1996). 
We then considered a cone of projected size 1.5/0.93 = 1.6 rioo- 
Adopting their typical k = 2, we then found that after a 2<x 
velocity cut, the fraction of particles with r > 1 .6 r2oo is now 
21.3%. This fraction is still marginally significantly larger than 
Biviano et al.'s fraction of 18% (assuming the same errors as 
above). We attribute this discrepancy to their variable k veloc- 
ity cut, which differs from our fixed k one. Wojtak et al. (2007) 
tried several interloper removal schemes and definitions (using a 
different ACDM cosmological simulation). Their local 3 cr cut 
leads to 20.4±1.7% of interlopers remaining within the virial 
cone. Given the quoted errors, the lower fraction of interlopers 
found by Wojtak et al. is marginally consistent with ours. 

This fraction of 23% of interlopers after the velocity cut is 
surprisingly close to the fraction of blue galaxies (i.e. galaxies 

19 The error is taken as their dispersion over the square root of their 
number of halos. 
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off the Red Sequence) observed within SDSS clusters, as Yang 
et al. (2008) find roughly 22% of blue galaxies within SDSS 
clusters of masses > 10 /2~'M . Admittedly, it is dangerous 
to match the dark matter distribution with the galaxy distribu- 
tion, since galaxies are biased tracers of the matter distribution. 
In fact, galaxies are biased relative to dark matter halos (e.g. 
Conroy et al. 2006), which in turn are biased relative to the 
dark matter particle distribution (e.g. Mo & White 1996; Catelan 
et al. 1998). If, in the end, the SDSS galaxies analyzed by Yang 
et al. are unbiased tracers of the dark matter distribution, then 
this close agreement would be expected if all blue galaxies are 
caused by projection effects. But if projections also pick up red 
galaxies in groups, then some blue galaxies would need to sur- 
vive within the virial sphere for the match to hold. However, the 
Yang et al. group finder is fairly efficient in separating groups 
along the line-of-sight, so we conclude that the fraction of blue 
galaxies within the virial sphere should be small. In other words, 
star formation appears to be strongly quenched when galaxies 
penetrate the virial spheres of clusters. 

When no velocity cut is performed, a maximum likelihood 
fit of the concentration of the projected NFW model to the pro- 
jected radii of a stacked cluster of nearly 300000 particles out 
to r2oo (/too) leads to a 14+7% (22+6%) underestimate of the 
true concentration parameter (Table 5, where most of the uncer- 
tainty comes from cosmic variance). Similar biases occur with 
the m-5 Einasto model. But after the velocity cut, these biases 
decrease by a factor two, and are no longer statistically signif- 
icant (Table 5). Moreover, the inclusion in such fits of a con- 
stant background as an extra parameter also strongly decreases 
the bias, even when the maximum projected radius is as low as 
r2oo (Table 5). In fact, inspection of Table 5 indicates that, for 
^max = '"too or 3 r2oo, the background (fixed or free) has a greater 
influence than the velocity filter in removing the bias on mea- 
sured concentration. Surprisingly, for R mdx = 3 rooo, a physically 
motivated fixed background added to the NFW model is slightly 
less effective in reducing the concentration bias than is a free 
background. 

When the maximum radius is 3 rioo and no velocity cut is 
performed, the NFW model with a free (respectively fixed) back- 
ground underestimates the concentration (Table 5) by 2+11% 
(15+9%). This insignificant (marginally significant) bias is 
caused by the strong decrease of the surface density profile once 
the Hubble flow is added to the peculiar velocities (Fig. 16). 
These small biases suggest that the fairly low concentration 
(C200 = 2.9+0.2) for the galaxy distribution in clusters found by 
Lin et al. (2004), who fit an NFW model with a free constant 
background to the distribution of projected radii in the range 
0.02 < R/r2oo < 2.5, but who did not make a velocity cut for lack 
of velocity data, is incompatible with true cluster concentrations 
of c = 4.0 at the 2 cr level. The lower concentration bias with 
the two-component model is expected, because the single com- 
ponent NFW or Einasto models cannot capture the flat surface 
density at large radii (Fig. 14), because other halos are projected 
along the line-of-sight. 

While a two-component model of halo (to infinity) + con- 
stant background is better able to recover the halo concentra- 
tion than a single-component model (Table 5), it is not wise to 
estimate the halo concentration from a two-component model 
with a halo term whose line-of-sight is limited to the sphere 
(Appendix B) plus a near constant background term arising from 
our universal interloper surface density model (eq. [ 1 2] or simply 
[13]): the single-component NFW captures better the total sur- 
face density profile than this halo+background model, especially 
if the maximum projected radius is beyond the virial radius, as 



the interloper surface density has a discontinuous slope at the 
virial radius (Fig. 14). On the other hand, the universal distribu- 
tion of interlopers in projected phase space might be useful to 
model the internal kinematics (hence total mass profile) of clus- 
ters of galaxies, where the full distribution of galaxies in pro- 
jected phase space is the sum of these interlopers and an NFW- 
like model projected onto the virial sphere. We are preparing 
tests of the mass/anisotropy modeling of clusters, groups, and 
galaxies (through their satellites) using this interloper model. 

We also performed 2D fits to individual halos of typically 
700 particles (not shown here). The dispersion of the concentra- 
tions were much larger (typically 0.16 dex) than the biases ob- 
tained from the stacked virial cone (typically 0.05 dex, i.e. 10% 
errors, see Table 5), which means that shot noise and cosmic 
variance dominate the bias caused by the Hubble flow. 

The line-of-sight velocity dispersion profile shows a concav- 
ity (in log-log) near the virial radius (Fig. 12), which is caused by 
the Hubble flow (Fig, 16). The velocity anisotropy profile recov- 
ered from this velocity dispersion profile, assuming the correct 
mass distribution, is close to the true anisotropy profile, with a 
slight, marginally significant, radial bias in the envelopes of clus- 
ters in comparison with the anisotropy profile recovered in 3D 
(Fig. 17), as was previously noted by Biviano (2007). 

In summary, the density profile of ACDM halos falls fast 
enough that the effects of the Hubble flow perturbing the stan- 
dard projection equations produce only small biases in compari- 
son with the shot noise of clusters with less than 1000 galaxies, 
as well as the large cosmic variance of the halos. 

These results have been obtained with the dark matter par- 
ticles of a cosmological /V-body simulation (with additional gas 
and galaxy components). They will need to be confirmed with 
future more realistic simulations of the galaxy distribution. 
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Appendix A: Projected mass, surface density and tangential shear of the Einasto model 

In this appendix, we derive an approximation to the surface density and projected mass (or, equivalently, projected number) profiles 
for the Einasto model. 

A. 1 . Projected mass profile 

For any density model, the projected mass is 
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(A.l) 



(A.2) 



where the second equality is obtained after reversing the order of integration. Equation (A. 1 ) is general, while equation (A.2) is only 
valid for models with finite total mass M«,. For the Einasto model of total mass M M , the 3D mass profile is 



M(r; m) = P 



3m, 2m | — 

r-2 



where P(a, x) = y{a, x)/F(a) is the regularized incomplete gamma function. The ratio 
M P (R; m) 



jj(R, m) = 



M(R; m) 



(A.3) 



(A.4) 



determined from equations (A.3) and (A.2), with equation (15), varies little, as seen in the right panel of Figure A.l. We fit again a 
two-dimensional fourth-order polynomial in m and u = log 10 (R/r^o) and find 



^a P x(", m) = dex (0.0001219 m 4 + 0.0007400 m 3 u - 0.003209 m 3 + 0.002976 mV - 0.01560 m 2 u 

+0.02966 m 2 + 0.0003307 mu z - 0.04434 mu 2 + 0.1273 mu- 0.1 149 m 
+0.001036 u A - 0.003133 u 3 + 0.1905 ir - 0.5241 u + 0.3525) . 



(A.5) 



(A.6) 



In the interval 3.5 <m< 6.5 and -2 < u < 2, equations (A.5) and (A.6) are accurate to better than 1.5% everywhere (0.23% rms). 
The projected mass of the Einasto model can thus be written 



M p (R;m) ^ yU apx (M, m) M(R; m) 
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(A.7) 
(A.8) 



where again u = log m (R/R-2) 



A.2. Surface density profile 

Inserting equation (15) into equation (1), the surface density of the Einasto model of total mass M and index m is 

UR;m) = — L\ — ;m\ , 

nr z _ 2 \r- 2 } 

3m " 1 xdx 



l (X ;m) = ^ exp(-2mjc 1/m ) 

y{3m,2m) J x v ' 



p(x; m) 



y{3m,2m) J x v ' Vx 2 - X 2 

Writing the dimensionless mass density as 
p(xr; m) (2m) 3m 



exp(-2mx 1/ffl ) 



(A.9) 
(A. 10) 



M(r_ 2 ;m)/(4jir 3 _ 2 ) my(3m,2m) 
where the second equality derives from equation (15), we can express the ratio of dimensionless surface to space densities as 



(A. 11) 
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(A.12) 
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Fig. A.l. Contours of log,,, 1 ?? (left, eq. [A. 12]) and ji (right, eq. [A.4], with eqs. [15], [A. 2] and [A.3]) for the Einasto model. 



where X = R/r-2- In the range 3.5 < m < 6.5 (spanned by ACDM halos in the redshift range < z < 3 according to Gao et al. 
2008) and -2 < log 10 X < 2, % varies little and regularly, as seen in the left panel of Figure A. 1 . We fit a two-dimensional 4th-order 
polynomial in m and u = log 10 (/?/r_2) to log 10 ??. We find 



<R(X,m) - K vx (u,m) 

Kpx(u,m) = dex (6.286 x 10~ 6 m 4 + 0.001 178 m\i - 0.0002251 m 3 + 0.001524 mV - 0.02427 m 2 u 

+0.0008538 m 2 + 0.001861 m u 3 - 0.02323 m u 2 + 0.1849 m u + 0.01577 m 
+0.0006014 m 4 - 0.01506 u 3 + 0.1056 u 2 + 0.3406 u - 0.2515) . 



(A.13) 



(A. 14) 



In the interval 3.5 < m < 6.5 and -2 < u < 2, equations (A.13) and (A. 14) are accurate to better than 0.8% everywhere (0.12% 
rms). 

The dimensionless surface density can then be written as 



T.(X; m) = 
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(A.15) 



or equivalently, with M (r_2) = P(3m, 2m) M^, where P(a, x) = y(a, x)/F(a) is the regularized incomplete gamma function and M a 
the total mass of the Einasto model: 
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(A. 16) 



Alternatively, the surface density profile can be, self-consistently, estimated from equation (A. 7) by differentiation over the 
projected mass profile, yielding after some algebra 
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(A. 17) 



where p is given in equation (A. 11). 

Equation (A. 17) has the advantage of providing an approximation for the surface density profile that is consistent with that of 
the projected mass profile. This is crucial for maximum likelihood estimation of concentration (and possibly Einasto index and 
background level). On the other hand, the accuracy of equation (A. 17) is about 5 times worse than that of equation (A.15). 
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A3. Tangential shear profile 

For any density model, the tangential shear measured by weak lensing can be written (e.g. Miralda-Escude 1991) 
X(R; m) - ~L(R; m) 

y t (R; m) = ; , (A. 1 8) 

■^crit 

where T.(R;m) = M p (R\m)l(jiR 2 ) is the mean surface density, while Z cr i t = c 2 /{4ttG) Ds/(Z)l£ ) ls) is the critical surface density, 
with c the velocity of light, and where D$, D^, and Z\s are the angular diameter distances between the observer and the source, 
the observer and the lens, and the lens and the source, respectively. Equation (A. 18) indicates that adding a constant term to the 
surface density (eq. [24]) has no effect on y t (this is the mass-sheet degeneracy). For the Einasto model, the tangential shear 
(eq. [A. 18]) is readily computed using equations (A. 15) with (A. 14) and (A. 8) with (A. 6). Figure A. 2 shows the subtle differences 
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Fig. A.2. Dimensionless tangential shear profile for the NFW model (black) and the m = 4 (red, long-dashed), 5 (green, short-dashed) and 6 (blue, 
dotted) Einasto models, using equation (A.18) with equations (A.9), (A.14), (A.15), (A.6), and (A.8). 

in the shear profile between the NFW and Einasto models of index m = 4, 5, and 6. While the tangential shear of the four models 
is indistinguishable in the wide range 0.8 < R/r^ < 10, there are potentially measurable differences at R > 10r_2 (at 100 r_2, 
the NFW shear is 1.5 times greater than that for the m-5 Einasto model) and possibly at R < 0.8 r_2 (as long as the weak linear 
approximation assumed for the measured shear to match the expression of y t of eq. [A.18] remains valid). 

Appendix B: Surface density and projected mass of the NFW model with lines of sight limited to a sphere 

In this appendix, we derive the surface density and projected mass (or, equivalently, projected number) profiles of the NFW model, 
with the lines of sight restricted to a sphere (which we conveniently choose as the virial sphere) instead of extending to infinity. 



B.1. Surface density profile 

In an analogous manner as the case with line-of-sight extending to infinity (eq. [1]), the surface density at projected radius R within 
the sphere of radius r max is 
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We now consider the case of the virial sphere: r max = r v . The surface density can then be written 



(B.l) 



(B.2) 
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where 

r p \x,c) = 



— f — 

'2-1 Jx (1 + 



2 1n2-l J x (1 + x ) 2 V* 2 -X 2 



1 



2 1n2-l 



(1 _X 2 ) 3 / 2 



cosh 



c + X 2 



(c+ 1)X 



1 Vc 2 - X 2 



(c+1) 1-X 2 



Vc^T(c + 2) | (-2c 3 -4c 2 -c + 2)(X-l) 
3(c+l) 2 SCc+l^V? 3 ! 



1 Vc 2 - X 2 1 



(c+1) X 2 - 1 (X 2 - l) 3 / 2 




cos 



c + X 2 



(c+ 1)X 



< X < 1 , 



X= 1 < c, 



1 < X < c , 



X = or X > c 



where equation (B.3) is found by inserting the NFW density profile (eq. [4]) into equation (B. 1). 
B.2. Projected mass profile 

For the NFW model, the projected mass within the virial sphere is 



where 



R r v 



2 n S S sph (S ; r v ) dS = M(r_ 2 ) M s J h | - 



r-2 r- 2 



M' p ph (X,c) = 2 



r x - 
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c)dY 



1 
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+ ln 
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cos 
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c + X 2 
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(B.3) 



(B.4) 



(B.5) 
(B.6) 

X = 0, 

< X < 1 and X < c , 

1 < X < c , (B.7) 

1 = X < c , 
X > c , 



where equation (B.7) was found by inserting equation (B.4) into equation (B.6). For X > c, one recovers the mass within the virial 
sphere. 

Appendix C: Maximum likelihood estimates 

In this appendix, we illustrate the maximum likelihood calculations that we have performed. 
Given parameters 6, and data points x the MLE is found by minimizing 



-\n£= -Y\np(xj\0), 

j 

where _£ = ^\pi. x jW) is the likelihood. 



(CI) 



C.1. Density profile 

The probability of measuring an object (galaxy or dark matter particle) at radius r in a spherical model of concentration c is 
47ir 2 \v{rf, c) + 



p(rj\c) = 



N(r m - dx ; c) - N(r min ; c) + Anb (4 ax - A n ) /3 ' 



(C2) 
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where v(R) and N(R) are respectively the density and number (proportional to mass) profiles, r m j n and r max are respectively the 
minimum and maximum radii, c is the concentration, while b is the constant density background. 

C.2. Surface density profile 

The probability of measuring a galaxy at projected radius R in a spherical model of concentration c and background b is 
2nRj h(Rj-,c) + -E bg ] 

P(Rj\c, b) = L - l — — , (C.3) 

Wnax; c) - N p (R mm ; c) + 7r£ bg (/?Lx - #y 

where 2,(R) and N P (R) are respectively the surface density and projected number (proportional to projected mass) profiles, R m j n and 
R m ax are respectively the minimum and maximum projected radii, c is the concentration, while Sbg is the constant surface density 
background. 

For the surface density profile S(7?) and the projected number (mass) profile N P {R), we use the formulae of Lokas & Mamon 
(2001) and of appendix A for the NFW and Einasto models, respectively. 

C.3. Distribution of interloper velocities 

According to equation (8), the distribution of interloper line-of-sight absolute velocities, vj = \v\ os J, is to first order the sum of a 
gaussian and a constant term: 

A exp \-v 2 l {la 2 )] + B 

p(vj\cr u A, B) = -— 1 I V 1,1 , (C.4) 

yfnJlA en erf ^/(cr; V2)J +kB 

where the denominator is found by ensuring C p(Vj) duj = 1 (eq. [12]), and where k is the maximum considered value of |ui os |/t>v (so 
k — 4 in Figs. 6). If A and B are expressed in virial units, then the denominator of equation (C.4) is the surface density of particles 
under consideration in virial units, which we directly measure from the simulation as 2 = (N/N v )/S , where N v is the number of 
particles within the virial sphere, while N is the number of particles in the radial bin (or within the full virial cone), and S is the 
surface of the radial bin (i.e. tt for the full virial cone). Hence, substituting for A — S (1 - kB')/[ yfnj2(Ji erf[£/(o"i V2)], we can 
write the probability of measuring an interloper absolute velocity as 

(I -kB') exp U 2 / (2^)1 
V^72cr ie rf [k/io-i V2)J 

where B' = fi/X. Then given the respective uncertainties e(cn) and e(B') in en and B' , we deduce the uncertainties in B and A as 

e{B) = ZeCB'), (C6) 

[/ dA \ 2 (dA\ 2 j2k 2 nor*S 2 e 2 (B') + (1 - kB') 2 { VSFtrifi - 2k exp [^/(2^f)lf eVi) 
<A) = J— ^Vi) + k^ e 2 (B') = X^ ■ l - , (C.7) 



[den) - ■ \dB'j ' v ' no-^8 2 
where S = erf ^/(cn V2)J. 



C.4. Practical considerations 

For one-parameter fits, we first search on a wide linear grid of equally-spaced 1 1 points for 0j, then we consider the three points 
with the lowest values of In £. and create a subgrid of 1 1 equally-spaced points (thus typically zooming in by a factor of 5), and 
iterate with finer subgrids until the two values of the parameter 6< with the highest likelihoods differ by less than 0.0001 or when 
the lowest - In £, decreases by less than 10~ 12 . We then obtain the 1 <x confidence interval fitting a cubic spline to the points below 
and above the best-fit parameter to solve for - ln_£ = - In -£ml + 0.5. 

For two-parameter (three-parameter) fits, we first search on a wide rectangular (cuboidal) grid of equally-spaced 1 1 points. Then 
we consider the rectangle (cuboid) obtained by searching for the lowest values of - In such that there are at least 3 different values 
for both (all three) parameters. We create a sub-grid in this rectangle (cuboid) with again 11x11(11x11x11) points, and iterate 
with finer subgrids until the pair of each of the two (three) parameters with the highest two likelihoods differ by less than 0.0001 or 
when the lowest - In X. decreases by less than 10~ 12 . We then obtain the 1 cr contour by considering those points in parameter space 
for which - lnX = - ln-CiviL + 1-15 (1.77), and then define as the minimum and maximum values for each parameter the extreme 
values in this contour. 
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